function [M, tList] = solveTransient(obj, tList, initM_vect)
    %solveTransient calculate M vector with given initial
    %conditions
    %   input:
    %       tList: double array, time array [s]
    %       initM_vect: [3 x 1] array, initial condition, M vector
    %   return:
    %       M: [3 x length(tList)] array, Mx = M(1, :); My = M(2, :); Mz = M(3, :);            
    %       tList: double row vector, time array [s]
    if nargin == 2
        initM_vect = [0;0;0];
    end
    p2 = inputParser;

    vectorChecker = @(x)validateattributes(x,{'numeric'},{'vector','real','nonnan','finite'});
    vector3Dchecker = @(x)validateattributes(x,{'numeric'},{'column','numel',3,'real','nonnan','finite'});
    p2.addRequired('tList',vectorChecker);
    p2.addRequired('initM_vect',vector3Dchecker);

    p2.parse(tList, initM_vect);
    
    if numel(tList) == 2
        tList = tList(1) : obj.parameters.maxStep/20 : tList(end);
    end
    
    options = odeset('MaxStep', obj.parameters.maxStep);
    BlochEquation = @(t, M) obj.blochEq(t, M, ...
                                        obj.parameters.T1, ...
                                        obj.parameters.T2, ...
                                        obj.parameters.OmegaX, ...
                                        obj.parameters.OmegaY, ...
                                        obj.parameters.OmegaZ, ...
                                        obj.parameters.Rp, ...
                                        obj.parameters.M0);
    sol = ode45(BlochEquation, tList, initM_vect, options);
%     [t,y] = ode45(BlochEquation, tList, initM_vect, options);

    % refine tList for better plot appearance
%     tList = refineGrid_1D(sol.x, 20);
    M = deval_fast(sol,tList);

    obj.solution.type = 'Transient';
    obj.solution.tList = tList;
    obj.solution.mVector = M;
    obj.solution.sol = sol;
end